Energetics and atomic mechanisms of dislocation nucleation in strained epitaxial layers 
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We study numerically the energetics and atomic mechanisms of misfit dislocation nucleation and 
stress relaxation in a two-dimensional atomistic model of strained epitaxial layers on a substrate 
with lattice misfit. Relaxation processes from coherent to incoherent states for different transition 
paths are studied using interatomic potentials of Lennard- Jones type and a systematic saddle point 
and transition path search method. The method is based on a combination of repulsive potential 
minimization and the Nudged Elastic Band method. For a final state with a single misfit dislocation, 
the minimum energy path and the corresponding activation barrier are obtained for different misfits 
and interatomic potentials. We find that the energy barrier decreases strongly with misfit. In 
contrast to continuous elastic theory, a strong tensile-compressive asymmetry is observed. This 
asymmetry can be understood as manifestation of asymmetry between repulsive and attractive 
branches of pair potential and it is found to depend sensitively on the form of the potential. 

PACS numbers: 68.55.Ac, 68.35.Gy, 68.90.+g 
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I. INTRODUCTION 

Emergence of misfit dislocations in heteroepitaxial sys- 
tems is a long-standing problem in the field of thin film 
growthiiSiMiSiSiLSiJ'i"'!!. Improving the physical proper- 
ties of semiconductor heterostructures requires control- 
ling the atomistic processes responsible for generation of 
defects. Thus, understanding the atomistic mechanisms 
of defect nucleation is crucially important for further 
progress in the field of heterostructure growth and struc- 
tural control of nanostructures. In addition, misfit dislo- 
cations represent an important problem in fundamental 
science. While a lot of information about the nature 
of dislocations has been obtained within the traditional 
continuum elastic theory, not near as much is known 
about the details of the underlying atomistic mechanisms 
through which dislocation nucleation occurs. 

Energy-balance arguments for the competition be- 
tween strain energy build-up and strain relief due to dis- 
location nucleation in mismatched epitaxial films lead to 
the concept of an equilibrium critical thickness. This is 
defined as the thickness at which the energy of the epitax- 
ial state equals that of a state containing a single misfit 
dislocation. It has been argued that dislocations should 
appear in the film when the thickness exceeds this critical 
valu o^i^i^ . The predicted critical value from this consider- 
ation, however, both from continuous elastic models'^ and 



from models incorporating layer discretenesses, is much 
smaller than the observed experimental value for the 
breakdown of the epitaxial state. This suggests that the 
defect-free (coherent) state above the equilibrium criti- 
cal thickness is metastableii and the rate of dislocation 
generation is controlled by kinetic considerations instead. 

The idea of strain relaxation as an activated process 
is supported by experimental results for the tempera- 
ture dependence of the critical thicknessii^iSii^ and it 
is the fundamental assumption in kinetic semi-empirical 
models^**. Physically, the lowest energy barrier for the 
nucleation of dislocations should correspond to a tran- 
sition path that initiates from the free surface (with or 
without defects). Such processes have been considered 
in a number of studies using continuum modelsi^iiSiii. 
However, it has been pointed out that surface steps and 
surface roughness that are not considered in the contin- 
uum models could play an important role for disloca- 
tion nucleationi2ii2i2£. Thus, atomistic studies are im- 
portant for a detailed understanding and determination 
of the possible mechanisms for defect nucleation in epi- 
taxial films. Although the importance of kinetic fac- 
tors in real experiments have already been emphasizedii 
and also investigated in numerical simulations of atom- 
istic models of the growth process^", direct determination 
of the transition path and corresponding energy barrier 
for misfit dislocation nucleation from an epitaxial film 



2 



has been much less explored, and they often require as- 
sumptions on the particular structure of the intermediate 
configurationSi. 

The actual stress relaxation processes starting from 
the epitaxial coherent state can occur along many dif- 
ferent transition paths. The path with the lowest acti- 
vation energy barrier at the saddle point corresponds to 
the true nucleation barrier for the generation of a mis- 
fit dislocation. For correct determination of this bar- 
rier, it is important to investigate different minimum 
energy paths (MEPs)^^, from the metastable coherent 
state to the incoherent state, without assuming a pri- 
ori any particular form of the intermediate configura- 
tions. We have recently carried out such a task which 
systematically explore the MEPs in the phase space of 
the system2iS4 based on a combination of the Repulsive 
Bias Potential (RBP)2^ and the Nudged Elastic Band 
(NEB) methods^. In the previous workS^, we consid- 
ered the case of a relatively large misfit of / = ±8 %. We 
showed that there is indeed a nonzero energy barrier for 
defect nucleation. Most importantly, however, we showed 
that both the mechanisms for the initiation of a misfit 
dislocation and the activation barrier exhibit a strong 
tensile-compressive asymmetry which is sensitive to the 
range of the interaction potential. A tensile-compressive 
asymmetry has also been found previousl}«22iSi in other 
contexts. 

In this work, we present a detailed systematic study of 
defect nucleation for the same 2D Lennard-Jones system 
as in Ref. ^23.. We consider strains in the range / — 
±(4 — 8) %, and intcrmolccular potentials with different 
ranges. 



II. MODEL 



equilibrium interatomic distance tq was set to a differ- 
ent value Tss, Tff and rfs for the substrate, film and film 
substrate interaction respectively. The parameter rg was 
varied to give a a misfit between lattice parameters as 

/ = (fff - rss)/rss- (3) 

For the film-substrate interaction, we set the equilibrium 
distance rfs as the average of the film and substrate lat- 
tice constants, i.e. rfs = (rff -|-rss)/2. A positive mis- 
match / > corresponds to compressive strain and neg- 
ative to tensile strain when the film is coherent with 
the substrate. Calculations were performed with peri- 
odic boundary conditions in the direction parallel to the 
film-substrate interface. For large systems, free boundary 
conditions gave qualitative similar results. In the calcu- 
lations, the two bottom layers of the five-layer substrate 
were held fixed to simulate a semi-infinite substrate while 
all other layers were free to move. Typically, in our cal- 
culations each layer contained 50 or more atoms. The 
central portion of the initial epitaxial film and substrate 
are shown in Fig.J^I- 

In the previous workS^ it was found that some fea- 
tures of dislocation nucleation are sensitive to the de- 
tailed form of the atomic potentials used. The results 
presented here are from systematic calculations for dif- 
ferent values of cut-off distances for the 5 — 8 potential 
(m = 8, n = 5). The advantage of this potential over the 
conventional 6 — 12 LJ potential is that it is intrinsically 
longer ranged. Thus, by imposing different cutoff radius 
Tc, one can study the influence of the range of the poten- 
tial on the nucleation of misflt dislocations. The other 
difference with respect to the 6 — 12 potential is a softer 
repulsive core. This will lead to a weaker anharmonicity 
and less asymmetry between the tensile and compressive 
strain situations. 



We consider a 2D model of the epitaxial film and sub- 
strate where the atomic layers are confined to a plane as 
illustrated in Fig. 1. Interactions between atoms in the 
system were modelled by a generalized Lennard-Jones 
(LJ) pair potentia!2^i that is modified to ensure that the 
potential and its first derivative go to zero at a predeter- 
mined cut-off distance ry. 



U{r) = V{r), r < tq; 



U{r) = V{r) 



ra 
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where 



V{r) = e 



n — m \ r 



ro 



n /ro 



, r>ro, 
(1) 

(2) 



r is the interatomic distance, e the dissociation energy 
and rg the equilibrium distance between the atoms. This 
potential for m = 12 and n = 6 is the same that has been 
used by Dong et al.^'^ in a recent simulation study. The 



III. METHOD 

The standard way of generating transition paths 
through Molecular Dynamics (MD) simulations^^ does 
not work well in cases where the probability for rare acti- 
vated events is small. There are now numerous methods 
which have been constructed to solve this fundamental 
problem. The MD technique itself has been augmented 
by various acceleration^^ and sampling schemes22i^. In 
addition, there is a class of methods that do not evaluate 
the dynamics directly but instead focus on a systematic 
search of transition paths and related saddle points for 
many-particle systems^it^^iSi^. 

We have recently introduced^A a particularly simple 
but efficient method called the Repulsive Bias Potential 
(RBP) method for transition path searching. In the RBP 
method, the potential energy of the system is augmented 
with a fixed, repulsive bias potential to make the ini- 
tial configuration unstable, but to keep the other nearby 
minima unaffected: 

C/tot(r,ro) = Uif) + Aexp{-[(r'- ro)/a?}- (4) 
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Here U{r) is the original potential energy surface of the 
system, which has been modified by an exponentially de- 
caying, spherically symmetric potential of strength A and 
range a which is centered at rp . When A and a have been 
chosen appropriately, forces computed from Eq. (1) can 
be used to displace the system from its initial state lo- 
cated at To to escape to a nearby minimum. This is done 
by applying total energy minimization to [/tot- 

With the RBP method implemented, the procedure of 
determining the transition path comprises several stages. 
First, the initial epitaxial state is prepared by minimizing 
the total energy of the system using MD cooling. In the 
MD cooling method, the energy is gradually minimized 
by setting the velocities v = whenever v and the force f 
on a particle satisfy the condition v-f < 0. Positions and 
velocities of the particles are obtained from numerical 
integration of the equations of motion using the standard 
leap-frog algorithm. Following this, the RBP is applied 
and the system is slightly displaced from the initial state 
(randomly or in a selective way to escape from harmonic 
basin) and then total energy minimization is applied to 
find a new minimum energy state. 

It is important to note that the RBP method can gen- 
erate many different final states depending on both the 
initial displacements and the exact form of the repulsive 
bias introduced. By making the repulsive bias sufficiently 
localized around the initial potential minimum, the final 
state energy depends only on the true potential of the sys- 
tem and not on the fictitious repulsive bias. In this work, 
we only consider the final configurations corresponding to 
the presence of a single misfit dislocation. Rather than 
trying random initial displacements, some knowledge of 
the dislocation generation mechanism is useful for expe- 
diting the process. 

We also find that the proper choice of initial displace- 
ments depends on the sign of the misfit. In the case of 
compressive strain, to get an ideal single dislocation lo- 
cated in the center of our sample, the optimal initial dis- 
placement corresponds to moving one atom in the middle 
of the first layer of the film from the film-substrate in- 
terface upwards by a small distance (0.04rss ). In case of 
tensile strain, the corresponding optimal initial displace- 
ment is a small displacement (0.04rss ) downwards for an 
atom located in the middle of the second layer in the film 
from the film-substrate interface layer. 

While the repulsive bias potential minimization can be 
used to generate the final state configuration containing a 
misfit dislocation, it does not yield the precise minimum 
energy path and the lowest activation barrier value for 
getting to this final state configuration. For this purpose, 
we use the Nudged Elastic Band (NEB) methodS^. This 
is an efficient method for finding the minimum energy 
path (MEP) , given the knowledge of both initial and 
final states. The MEP is found by constructing an initial 
set of configurations (images) of the system between the 
initial and final states. This set is then allowed to relax 
to the true set representing the MEP. 

An initial guess for the images in the NEB is usually 



obtained by interpolating the particle configurations be- 
tween the final and initial states. For the present applica- 
tion, however, we find that this often leads to numerical 
instabilities due to the strong hard core repulsion of the 
LJ potentials and fail to converge to the true MEP. To 
circumvent this problem, we use the set of configurations 
generated in moving to the final state in the presence of 
the repulsive bias as the initial input in the NEB. This 
leads to fast convergence in the NEB method without 
the instabilities encountered in the linear interpolation 
scheme. 



IV. RESULTS 

For epitaxial films above the equilibrium critical thick- 
ness, the relaxed state with a nonzero density of misfit 
dislocations which partially relieves the strain energy in 
the film is expected to have a lower energy. However, if 
this configuration is separated from the coherent state by 
a finite energy barrier AE, the film will remain coherent 
unless defects are nucleated allowing to overcome this en- 
ergy barrier. This barrier could be finite even when the 
relaxed state has already a lower energy than the epi- 
taxial state. Thus the experimentally observed critical 
thickness can be much larger than the equilibrium value 
depending on the kinetics of defect nucleation. Our pre- 
liminary results^-2i^ showed a large variety of relaxation 
processes, including single dislocation nucleation, mul- 
tiple dislocations, dislocations with different core struc- 
tures, and dislocations nucleating on different depth in 
the film, which can be characterized by their different 
activation energies and energies of the final incoherent 
states. In this work, we focus on the nucleation and 
MEP leading to a final state containing only a single mis- 
fit dislocation with core located near the film-substrate 
interface. To simplify the discussions, we will present in 
this section only the results for the 5 — 8 potential with a 
cutoff radius of Tc = l-5rss, and lateral size L = 50, corre- 
sponding to 50 atoms per layer. These results allow us to 
arrive at a simple physical picture for the nucleation pro- 
cess of the misfit dislocation. The results with different 
parameters for the intermolecular potential and different 
size of the system are qualitatively similar. They will be 
presented in a later section. 



A. Mechanisms of relaxation 

Relaxation of strain with dislocation nucleation is a 
complex process involving motion of many particles in- 
side the system. The transition from coherent to dislo- 
cated state considered in this paper is analogous to strain 
relaxation in a real heteroepitaxial sample under anneal- 
ing conditions. Experiments show that heating is a es- 
sential prerequisite for such relaxation to occuriii*!^. This 
fact shows that nucleation of dislocation represents typ- 
ical activated process with a nonzero activation barrier. 
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Our calculations with NEB confirm this conclusiori2i2a. 
For both the compressive and tensile strain cases, we find 
the presence of a finite activation barrier AE along the 
MEP leading from the initial epitaxial state to the final 
state with a single misfit dislocation in the film substrate 
interface. To allow for comparison of different cases and 
extraction of the basic physics involved, we introduce the 
definition of the reaction coordinate S. This is defined 
as the accumulated displacement of the system along the 
MEP in the multidimensional configuration space. Math- 
ematically, the reaction path coordinate for a given con- 
figuration (image) along the MEP is defined as 



M N 

^*^ = E\E(^™-C"')V^ (5) 

m— 1 1 i—1 

Here M is the label for the configuration (image ) under 
consideration, and i is the index for the different parti- 
cles in the system (i=l to N). In Fig.© and Fig.Q, we 
show typical snapshots of configurations along the corre- 
sponding MEP for compressive and tensile strain cases 
respectively. In all cases the initial state was an epitaxial 
film with a coherent interface and the final state con- 
tained a single dislocation with its core located in the 
interface layer. The final state is characterized by the 
presence of an adatom island on the surface of the film 
in the case of compressive strain and a vacancy island in 
the tensile case. The number of adatoms (or vacancies) 
in the island exactly corresponds to the number of layers 
in the film. Such form of the final state is determined by 
the geometry of the misfit dislocation, as the one extra 
atom is added or removed from (or inside) each layer to 
relax the strain. 

An important property of the NEB method is that it 
usually converges to the MEP nearest to the initial trial 
trajectory. Thus by changing the initial input path, we 
were able to investigate several different mechanisms of 
relaxation^^'^'*. These mechanisms differ from each other 
mainly by the level of coUectiveness in the displacement 
of the particles from the coherent state position. For each 
given set of parameters, we identify the lowest activation 
barrier. The particular kind of mechanism leading to the 
lowest activation barrier depends on the parameters of 
model (misfit, cut-off radius of potential etc.). We find 
that for all the systems that we have studied, the mech- 
anisms leading to the lowest activation barrier belong to 
one of the two categories described below. 

The first mechanism describing the transition from the 
initial coherent state to the final state with a misfit dis- 
location at the film substrate interface corresponds to a 
successive sliding along the edges of a triangle. The sad- 
dle point configurations corresponding to this mechanism 
for the tensile and compressive strain cases are shown in 
Fig. 10^) and Fig. Qd) respectively. We see that in this 
case the displacements of the atoms have a collective be- 
havior, with two edges of a triangle successively sliding 
up or down (one by one). Eventually, an adatom island 



or a vacancy island is created on the surface of the film. 
The highest saddle point can correspond either to the 
sliding of the first or the second edge. We refer to this 
as the glide mechanism since the motion of the disloca- 
tion after it is nucleated follows the path referred in the 
literature as dislocation glided. For the tensile strained 
film, the glide mechanism always yield the lowest activa- 
tion barrier. While for the compressively strained film, 
the mechanism leading to the the lowest activation bar- 
rier depends actually on the magnitude of the misfit. For 
small misfit ( / ^ 8%), the glide mechanism is again 
the one leading to the lowest activation barrier. This is 
drastically different from the climb mechanism reported 
earlier— for a misfit of 8% in a compressively strained 
film. 

The second mechanism correspond to successive relax- 
ation of layers. This is the preferred mechanism for a 
compressively strained film with large misfits (/ ^ 8%). 
The saddle point configuration corresponding to this 
mechanism for the compressive strain of 8% misfit is 
shown in Fig. 10];). In this case, the core of the disloca- 
tion first appears at either the second or the third layer of 
the film and then successively moves down from layer to 
layer to the film-substrate interface. The displacement of 
the particles have a very localized character in this kind 
of mechanism. We refer to this as the climb mechanism 
since the motion of the dislocation after it is first nucle- 
ated in this case corresponds to what is known in the 
literature as dislocation climb^. For intermediate values 
of compressive strain, the situation is more complicated, 
as the two mechanisms are competitive in energy costs. 
The actual MEP in this case is better described by a 
mixture of the climb and glide mechanisms. 



B. Activation energy of dislocation nucleation 

The most important characteristic of a particular re- 
laxation process through nucleation of a misfit disloca- 
tion is its activation energy AE. The activation barrier 
is calculated as the difference between the total energy 
of the initial state and that of the saddle point config- 
uration. As can be seen in Fig. ||2Jl, corresponding to 
the compressive strain case, there may exist many saddle 
points along a given MEP . The activation barrier is de- 
termined by the highest energy saddle point. The results 
for AE vs. the number of layers in the film are presented 
in Fig.©. 

For the tensile strain case, we find that the process 
leading to the nucleation of misfit dislocation and sub- 
sequent motion along the MEP is always through the 
glide mechanism. The activation barrier decreases with 
increasing magnitude of misfit. Also, at large misfits, 
the activation barrier decreases significantly as the film 
thickness increases, leading to an essentially negligible 
activation barrier. This was verified directly through in- 
dependent MD simulation at finite temperatures where 
the misfit dislocation is easily generated spontaneously. 
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For the compressive strain case, except at 4% mis- 
fit and small thickness (less than six layers), the barri- 
ers are higher than the corresponding tensile strain case 
with the same magnitude of misfit. Again, there is a 
strong decrease in /S.E with increasing magnitude of mis- 
fit. In contrast to the tensile strain case, the activation 
barrier tends to level off with increasing film thickness. 
The other striking difference from the tensile strain case 
is that the mechanism corresponding to the movement 
along the MEP in this case can either be the glide mech- 
anism as in the tensile strain case, or the qualitatively 
totally different climb mechanism involving layer by layer 
distortion as discussed in the last section. This new climb 
mechanism occurs for large misfits (/ ^ 8%). 



V. SIMPLE PHYSICAL PICTURE FOR THE 
NUCLEATION PROCESS 

As shown in the last section, the mechanism leading 
to the nucleation of a misfit dislocation starting from 
the epitaxial coherent state and the subsequent motion 
along the MEP to the final state is fairly complicated, 
and depends sensitively on the sign and magnitude of the 
misfit (tensile or compressive strain) , and thickness of the 
film. With this rich set of data, it is important to have 
some simple qualitative understanding of the results. 

First of all, it is easy to understand the origin of the dif- 
ference between the tensile and compressive strain cases. 
In a harmonic elasticity theory, the activation barrier 
would depend only on the magnitude and not the sign 
of the strain. The tensile-compressive asymmetry thus 
originates from the strong anharmonicity of the interac- 
tion potential, particularly in the steeply rising repulsive 
core. This is confirmed by our results shown in Fig.Q 
showing that the difference of IS.E for the tensile and com- 
pressive cases grows monotonically as the misfit increases 
in magnitude. This is also confirmed in our similar stud- 
ies using the conventional 6 — 12 LJ potential as shown 
in Fig.®. Since the 6 — 12 potential is considerably 
steeper in the core region, the anharmonicity is stronger 
and the resulting tensile-compressive asymmetry is even 
more pronounced. 

The other general trend is the strong decrease of the 
activation barrier with increasing misfit. This is true for 
both the tensile and compressive cases (Fig. O. It re- 
mains true even when the mechanism leading to the nu- 
cleation has changed character from a glide nature to a 
climb nature as in the case of large compressive strain. In 
our previous work^?:, we have analyzed the contribution 
to the activation barrier from the intralayer and inter- 
layer bond distributions at the saddle point. Here we 
will introduce the same physical arguments in terms of 
the conceptually simpler quantity of reaction coordinate 
defined earlier in Eq. jSJ. Let S represent the dimension- 
less reaction coordinate along the MEP leading from the 
initial coherent state through the saddle point to the fi- 
nal state containing the misfit dislocation. For the initial 



stages of small displacement with 5' -C 1, the simplest 
leading representation of the MEP can be expressed in 
the form 

E=''-S^~ls\ (6) 

In the equation above, the first term gives the energy 
rise towards the saddle point from the initial displace- 
ments from the coherent state necessary to nucleate the 
dislocation. It originates mainly from the stressing of the 
interlayer bonds which are fully relaxed in the initial co- 
herent epitaxial state. Because of this initial relaxation, 
there is relatively little dependence of the coefficient a 
on the misfit. The second term represents the release 
of the intralayer strain energy from the displacements of 
the atoms. Clearly, the coefficient h is strongly dependent 
on the magnitude of the misfit. Whether it is tensile or 
compressive, the higher the magnitude of the strain, the 
larger is the lowering of the strain energy. Hence the co- 
efficient b should be a monotonically increasing function 
of the magnitude of the misfit. It follows simply from 
Eq. © that the activation barrier Ai? is given by the 
expression 



Thus, the activation barrier always decreases with in- 
creasing magnitude of the strain, whatever the actual 
initial strain release mechanism and nature of the sad- 
dle point configuration. Furthermore, the expression in 
Eq. © predicts that the saddle point should occur at 
the reaction coordinate 5o = a/h which again decreases 
monotonically as the misfit magnitude increases. This is 
supported by our results as shown in Fig.(|7|l. 

In general, the initial cost of energy in creating the dis- 
tortion for the dislocation in the glide mechanism is lower 
for a tensile than compressive strain. This is due to the 
fact that for the compressive strained film, the initial dis- 
tortion required for creating the dislocation core always 
involve breaking of bonds to lower the coordination num- 
ber. On the other hand, for the tensile strained film, no 
breaking of bonds is necessary in the glide mechanism for 
the nucleation of the dislocation. Thus the glide mecha- 
nism is always preferred for the tensile strained film. For 
the large compressive strain, the energy cost involved in 
nucleating a dislocation core is comparable for the glide 
and climb mechanism, and the two processes are compet- 
itive. 

The dependence of the activation barrier on the film 
thickness is more complicated and is rather different for 
the tensile and compressive strains. For the large com- 
pressive strain case where the MEP corresponds to the 
climb mechanism, the behavior is fairly easy to under- 
stand as the saddle point involves a rather localized dis- 
location in the surface layers, so obviously the activation 
barrier would have a very weak dependence on the film 
thickness as observed in our numerical study. For the 
glide mechanism, both the initial rise in energy and the 
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release of the strain energy leading to the saddle point 
configuration are dependent on the film thickness and ac- 
cording to Eq. ((nj , it is hard to predict any universal de- 
pendence of the activation barrier on the film thickness. 
Indeed, both a levelling off (for compressive strain) and 
a strong decrease in the activation barrier as a function 
of the film thickness have been observed. 



VI. SIZE AND POTENTIAL DEPENDENCE 

The results presented in the previous sections are all 
for a 5 — 8 short ranged L-J potential with a cutoff set at 
l.Srss- The size of the system was set at L = 50 particles 
per layer. We have also performed similar calculations 
for different set of parameters in the potential as well 
as for different sizes to investigate the size and potential 
dependence of our results. We find that the results with 
different interatomic potentials and sizes of the system 
arc qualitatively similar, although differing in details. We 
present some of these results in this section. 

In Fig. ((HJ, the activation energy barrier is plotted 
against the film thickness for a system size of L = 20 
and a short ranged potential as in the previous sections 
for two values of the magnitudes of misfit at |/| = 5% 
and 8%. The results are very similar to that presented in 
Fig. 10. The only limitation for the smaller sample size 
is that one cannot accurately study the cases of smaller 
misfit as the addition or removal of a single atom from a 
layer would overshoot the strain release mechanism. 

In Fig. (|5J), we show the results of activation energy 
barrier vs film thickness for system size L = 50 and a 
5 — 8 LJ potential as before but this time with a longer 
range with cutoff set at 2.1rss. Again, the results are 
qualitatively similar to that presented in Fig. The 
tensile and compressive asymmetry is stronger for this 
longer ranged potential, particularly at the smaller misfit 
values. This could be also related to the stronger size 



effects for the longer ranged potential. 



VII. CONCLUSIONS 

We have developed a general scheme of identifying min- 
imal energy paths for spontaneous generation of misfit 
dislocation in an epitaxial film and studied the energet- 
ics and atomic mechanisms of stress relaxation using a 
two-dimensional model. This approach requires no a pri- 
ori assumptions about the nature of the transition path 
or the final states. A nonzero activation barrier for dis- 
location nucleation is found in the minimum energy path 
from coherent to incoherent state. We find that the en- 
ergy barrier decreases strongly with misfit. The nucle- 
ation mechanism from a flat surface depends crucially 
on whether we start from a tensile or compressive ini- 
tial state of the film. This asymmetry originates from 
the anharmonicity of the interaction potentials which 
leads to qualitatively different transition paths for the 
two types of strains. The present method can also be 
extended to three-dimensional models with more realis- 
tic interaction potentials. Preliminary calculations for a 
three-dimensional Lennard- Jones system and the Pd/Cu 
and Cu/Pd systemg-^^ with the Embedded Atom Model 
potentials'^^ confirms the effectiveness of the method in 
three dimensions. 
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FIG. 2: Minimal Energy path for compressive strain / = 
+8% as a plot of energy barrier AE vs reaction coordinate S. 
Snapshots configurations (a), (b) and (c) correspond to the 
labels in the energy profile (top right). Closed line in (c) is 
the Burgers circuit around the dislocation core. The energy 
barrier is in units of interatomic potential strength e and the 
reaction coordinate S is in units of equilibrium distance rss- 
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FIG. 1: A 2D model of the epitaxial film and substrate show- 
ing the particle configurations in the coherent state. The two 

layers at bottom arc held fixed while all other are free to move. 
Filled circles represent the the epitaxial film and open circles 
the substrate. 



• •••••••• -2 

>••••••••• QJ 

OOOOOOOOO <1 
)OOOOOOOOC -4 

OOOOOOOOO 
)OOOOOOOOC 

0.0 0.1 0.2 

S 



•••••••b 

••• •••••• 

•••••• 

• ••• ••••• 

••••• 

OOOOOOOOO 
)OOOOOOOOC 
OOOOOOOOO 
)OOOOOOOOC 



• • • • ( c 
»••••••••• 

• • • • 

• • • • 4 

O Of) O O JZ O O o 
) o o o O (» » o o c 

OOOOOOOOO 
)OOOOOOOOC 



FIG. 3: Minimum Energy path for tensile strain / = —8% as 
a plot of energy barrier AE vs reaction coordinate S. Snap- 
shots configurations (a), (b) and (c) correspond to the labels 
in the energy profile (top right). Closed line in (c) is the Burg- 
ers circuit around the dislocation core. The energy barrier is 
in units of interatomic strength e and the reaction coordinate 
S in units of equilibrium distance rgs- 
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FIG. 4: Saddle point configurations for different mechanisms 
of stress relaxation (a) glide mechanism for tensile strain (b) 
glide mechanism for compressive strain (c) climb mechanism 
for compressive strain. Filled circles represent the the epitax- 
ial film and open circles the substrate. 




FIG. 5: Energy barrier A_E (in units of e) as a function of 
film thickness (number of layers) for different misfit values. 
Squares symbols correspond to / = ±4%, stars to f — ±5%, 
triangles to f — ±6%, and circles to / = ±8%. Solid and 
dotted lines correspond to compressive / > and tensile / < 
strains, respectively. 
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FIG. 6: Energy barrier AE (in units of e) as a function of 
film thickness (number of layers) at misfit 5%, for the 5 — 
8 (squares) potential and 6 — 12 (circles) potential (cut off 
1.5rss). Solid and dotted lines correspond to compressive / > 
and tensile / < strains, respectively. Here the system size 
is L = 20. 



-10 




FIG. 7: Energy profile of the minimum energy path for (a) 
compressive and (b) tensile strain and for different misfits. 
Energy in units of e and S in units of equihbrium distance 
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FIG. 8: Energy barrier AE (in units of e) as a function 
of film thickness (number of layers) for smaller sample size 
(20 atoms per layer) and different misfit values for the 5 — 8 
potential: / = ±5% ( stars), and / — ±8% (circles). Solid 
and dotted lines correspond to compressive / > and tensile 
/ < strains, respectively. 
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FIG. 9: Energy barrier AE (in units of e) as a function of film 
thickness (number of layers) for different misfit values for long 
ranged 5 — 8 potential (cut off 2.1rss): / ~ ±4%(squares), / = 
±5% (stars), / = ±6% (triangles), and / = ±8% (circles). 
Solid and dotted lines correspond to compressive / > and 
tensile / < strains, respectively. 



